文章同步發表至 Medium

前兩篇介紹了如何使用 NetTopologySuite 來讀寫 Shapefile 檔案，這篇要介紹的是最後一種座標轉換的方法——利用 NetTopologySuite 讀取 WKT 之後，搭配 ProjNet 進行座標轉換。

前置作業

還記得前幾篇所說的 WKT 以及座標系統嗎? 為了讓程式知道我們要轉換的是哪一個座標系統，我們要先取得說明資訊的 WKT，可以從 epsg 提供的網站找到，我自己都是選擇第二個選項的 WKT 複製：

利用 ProjNet 進行座標轉換

以下是 NetTopologySuite 2.0.0 以後的寫法，我自己在寫的時候因為使用了更早之前的方法所以一值遇到問題。詳細資訊可以參考這篇討論：https://gis.stackexchange.com/a/294263

// 座標系統的 WKT
const string Wgs84Wkt = "GEOGCS[\"WGS84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
const string Twd97Wkt = "PROJCS[\"TWD97 / TM2 zone 121\",GEOGCS[\"TWD97\",DATUM[\"Taiwan_Datum_1997\",SPHEROID[\"GRS 1980\",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"3824\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",121],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",250000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"3826\"]]";

Console.WriteLine("請輸入經度：");
var lon = Console.ReadLine();

Console.WriteLine("請輸入緯度：");
var lat = Console.ReadLine();

var wktBeforeConvert = $"POINT({lon} {lat})";
Console.WriteLine($"轉換前的 WKT 為：{wktBeforeConvert}");

// 將 WKT 讀成 Geometry
var wktReader = new WKTReader();
var originPoint = wktReader.Read(wktBeforeConvert);

var coordinateSystemFactory = new CoordinateSystemFactory();
var coordinateTransformationFactory =  new CoordinateTransformationFactory();
    
var from = coordinateSystemFactory.CreateFromWkt(Wgs84Wkt);
var to = coordinateSystemFactory.CreateFromWkt(Twd97Wkt);

var fromCoordinateSystems = coordinateTransformationFactory.CreateFromCoordinateSystems(from, to);
var convertedPoint = Transform(originPoint, fromCoordinateSystems.MathTransform);

var wktAfterConvert = convertedPoint.AsText();
Console.WriteLine($"轉換後的 WKT 為：{wktAfterConvert}");

static Geometry Transform(Geometry geom, MathTransform transform)
{
    geom = geom.Copy();
    geom.Apply(new MTF(transform));
    return geom;
}

sealed class MTF : ICoordinateSequenceFilter
{
    private readonly MathTransform _mathTransform;

    public MTF(MathTransform mathTransform) => _mathTransform = mathTransform;

    public bool Done => false;
    public bool GeometryChanged => true;
    public void Filter(CoordinateSequence seq, int i)
    {
        // 如果有需要轉換 z 座標的話可以把註解打開
        
        double x = seq.GetX(i);
        double y = seq.GetY(i);
        // double z = seq.GetZ(i);
        
        _mathTransform.Transform(ref x, ref y);
        // _mathTransform.Transform(ref x, ref y, ref z);
        
        seq.SetX(i, x);
        seq.SetY(i, y);
        // seq.SetZ(i, z);
    }
}

輸出的結果如下：

可以利用線上的轉換工具幫助我們檢查輸出的結果是否正確：

或是直接利用 PostgreSQL 轉換看看：

References

• StackExchang - How to perform transformations using projnet
• StackExchang - Transforming point using nettopologysuite
• epsg.io